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The 2nd order Runge-Kutta scheme Langevin simulation of unquenched QCD in pseudofermion method 
derived from our general theory shows a behaviour as a function of the Langevin step t better than the 
Fukugita,Oyanagi,Ukawa's scheme. 



1. Introduction 

In order to incorporate fermions in the lattice 
QCD, Hybrid Monte Carlo(HMC) is the stan- 
dard method 0], but in the case of large lattice, 
the method is not free from problems caused by 
rounding errors^]. The Langevin simulation is 
another possible method to simulate fields includ- 
ing fermions with or without gauge fixing. 

The discretized Langevin simulation results de- 
pend on the step size t, and the incorporation of 
higher order corrections in the step function to 
reproduce the Fokker-Planck distribution is the 
problem. Drummond et al|3| showed a method to 
derive higher order corrections up to second or- 
der in the step in U(l) spin model and in SU(N) 
lattice gauge theory. We obtained the Langevin 
step function up to third order in flat space and 
up to second order in curved space 0]. 

A Langevin simulation using the Runge-Kutta 
type algorithm including dynamical quark loops 
and corrections up to second order in t was 
proposed by Batrouni et alQ and, Ukawa and 
Fukugita§. 

An extensive full-QCD Langevin simulation 
was performed by Fukugita, Oyanagi and 
Ukawa fit] in partial 2nd order bilinear noise 
scheme. They reported that the 2nd order 
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pseudo-fermion scheme is worse than the partial 
2nd order bilinear noise scheme in a simulation of 
4 3 x 8 lattice. 

In this work, we extend the new algorithm 
of Langevin simulation to an algorithm in the 
canonical coordinate space of Lie algebra of 
SU(N). We then apply the method to SU(N) lat- 
tice gauge theory including pseudofermion and 
discuss possible improvements. 

2. Random walks in the space of Lie groups 

In the space of SU(N) Lie groups, a group ele- 
ment g is expressed as g = exp(x • A) where Xi's 
are the canonical coordinates of the first kind, and 
Aj's are antihermitean matrices which span space 
of su(N) Lie algebra as [Aj, Aj] = CijkXk and are 
normalized as Tr(\\\j) — Sij . 

The right differentiation V is defined for any 
f{g) as e x v f(g) = f{ge x x ) , where Vj's follow 
the same commutation relation as Aj's and 

e -X.V 6 inv {g _ go) = ginv (g _ ^X-X) (1) 

The Fokker-Planck equation in the space of Lie 
group is d t (j> = Kcj) , where K = V fc (V fc - u k ) , 
and <j) is normalized with the Haar measure. 

If Ui — ViS, then <p tends to its asymptotic 
distribution const- e s . The formal solution </> with 
the initial condition <fi(g) = 5 mv (g — go) is 

4> = e tK 5 inv (g-g ). (2) 
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The Langevin simulation algorithm for the 
random walk in the Lie algebra space is defined 
as 



„tK jrinv 



6 mv (g~ga) = (6 mv (g-goe x - A )), 



(3) 



where X is a function determined from s = 
\/t, and all values of Uj together with their 
higher derivatives evaluated at <7oj and indepen- 
dent Gaussian random variates £'s with variance 
2 and mean 0, and (...) denotes an averaging over 
the Gaussian distribution. 

Thus the problem is to find the Xi such that 



(g\eMtK)\g ) = (g\(e 



\9o) 



(4) 



where K = V-(V — u(g)), and note that operators 
are normal ordered, i.e., all V's are set in the 
leftmost position relative to all g's on the l.h.s. of 
the equation. 

In order to get X, let Xi be expanded as 



(5) 



n=0 k=0 



where /ifc(£) is a generalized k-th order Hermite 
polynomial. First, put exp(tK) in normal order 
as mentioned above. Then starting from Xi = 
s£i + • • •, and considering operator equations for 
each order 0(s n ), and then for given n, 0(V m ) 
in larger-m-first order, we proceed perturbatively 
and obtain f( n ,k)'s- We can prove this procedure 
works satisfactorily, i.e., a newly fixed ft n ,k) does 
not affect lower order equations. 
Our 0(s 6 ) algorithm is as follows, 
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+12£jVfcUj VfcWi + 48£jVjttfcVfcUi 
+48£j<UfeVfcVj'Uj + 
+96^V :( V 2 M l + 24cy fc cy m ^V fc u m 
+2Actjk£,iUjViUk + 48cijk£iVjViUk 
~2Acjki^jUkViUi + 4&Cjki€jVkViiH - N 2 £i) 



-—(-ZONujVjUi - 10VV-V 
60 



-lOujVjUfcVfcUi - 30NV 2 m 
— 20VjWfeVjVfcUi — lOujUfcVfcVjiii 
-20n i V i V 2 u i - 10V 2 V 2 it l - 10ci jk VjUiViUk 
-5ci 3k u.jV 2 Uk - ^CijkUjUiViUk 

'^CijkCljrn'U'l^ m^k ^^ijk^lkm^l^ j^m 

-54N 2 u t ) (6) 

In the case of SU(3), exponentiation of the ma- 
trix x ■ A is done analytically as follows. 



e lZl P 1 +e lZ2 P 2 +e lZ3 P 3 . 



(7) 



with the eigenvalue given by the Cardano 

method, and Pi as the projection operator on its 
eigenspace. The unitarity of the matrix is pre- 
served in very high accuracy. 

2.1. QCD with pseudofermion 

In order to incorporate quark loops in the vac- 
uum, we include in addition to the gauge field U, 
a pseudofermion field <f> which is a complex-scalar, 
and define the partition Z = J dcf>*d(j)dUe" s . 

The action is 



S — 



Sfe 



where 

Sgauge = "T7 ^ Re(TrU p ) 



(8) 



(9) 



The matrix elements for Wilson fermions are 
calculated by using hermitian 7 matrices and the 
hopping parameter k. 

After the preconditionings the fcrmionic action 
can be replaced by the even-site fcrmionic ac- 
tionllOtl 



Sfermi- € 0e 



(10) 



J "1 



where M = 1 - n 2 DD and D xy = J2^{ 
(1 - ^)U X ^S 

The 0(s 4 ) Runge-Kutta algorithm X for the 
SU(N) link variables U is 

X = ^{-st, + « 2 2w) + s£ + \s 2 {u + u) (11) 
where u = u(ge^~^ s u ^' x ). The velocity field is 

^ijj, — ^ gauge H - ^"/i^/ermi- (^^) 
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where, is the right differentiation with respect 
to the i-th canonical coordinate of U Xtli . 

The function X for the pseudofermion is sim- 
ilar but the N dependent term is absent, u$ = 



AfMt 



and 



-2(f)* 



MMt 



3. Numerical Results 

In the second order Taylor scheme Langevin 
simulation of SU(2) on 8 4 lattice, we found the de- 
pendence on t of the data simulated in the space 
of Lie algebra is much better than that in the 
space of S 3 §. 

Simulations of SU(2) single degree of freedom, 
and of 4 4 and 8 4 lattices were performed in the 
Taylor scheme and in the Runge-Kutta scheme. 
In the single degree of freedom both showed sim- 
ilar dependence, however in the 4 4 lattice, the 
behaviour of the Taylor scheme was worse than 
the Runge-Kutta scheme at large t||. 

In the case of SU(2) single degree of freedom, 
the third order Taylor scheme did not improve the 
second order scheme, but Runge-Kutta scheme 
that incorporates the third order term is expected 
to improve the dependence on t in the lattice sim- 
ulation. 

We checked the results of SU(3) by compar- 
ing with the Monte Carlo simulation, in the 
Cabibbo-Marinari 's scheme. The dependence 
of (Up) = (Re(\TrU p )) on k in the 4 4 lattice 
in the 2nd order Langevin simulation of full- 
QCD with pseudofermion is shown in Fig.l. At 
t = 0.05, the agreement to the HMC result is bet- 
ter than that of Fukugita, Oyanagi and UkawajQ 
and at t = 0.01 they are consistent within statis- 
tical errors. 

4. Outlook 

General algorithm for obtaining higher order 
weak Taylor approximations to the Langevin step 
algorithms with respect to t, improved version in 
practical use for a finite fixed order was formu- 
lated, (cf. Nakajima || for full order weak Taylor 
algorithm in recursion formula) 

We try to obtain a local 0(t 3 ) Runge-Kutta like 
algorithm in SU (3) and aim to use it for investi- 
gation of its performance in full QCD simulation. 



0.6 



< Up > 4 4 = 5 



0.5 - 



0.4 



0.3 



1 


1 1 

t=0.01 


l 




ref.7 t=0.01 


A 




t=0.05 


— i ~ 




ref.7 t=0.05 


X 


r 


} *iff 




1 


X 

i i 


i 



0.05 



0.1 0.15 

K 



0.2 



Figure 1. Dependence of Wilson loop Up = 
Re(^TrUp) on the hopping parameter k. The 
triangle at k = 0.15 is the result of HMC[1]. 
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